package com.itfreer.gis.slice;

import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconst;
import org.gdal.gdalconst.gdalconstConstants;
import org.gdal.osr.SpatialReference;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;

import java.io.File;

public class ImageSlice {

    public static void main(String[] args) {
        gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
        gdal.AllRegister();
        tifSliceTest("D:\\shp\\tif\\LGQ_18.tif", "D:/shp/tif/slice", 14, 14);
    }

    public static void tifSliceTest(String fileName, String outputDir, int minZoom, int maxZoom) {

        // 获取影像对象
        Dataset dataset = gdal.Open(fileName, gdalconstConstants.GA_ReadOnly);
        if (dataset == null) {
            System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
            System.err.println(gdal.GetLastErrorMsg());
            return;
        }

        double[] ori_transform = dataset.GetGeoTransform();

        int xCount = dataset.getRasterXSize(); // 栅格尺寸宽
        int yCount = dataset.getRasterYSize(); // 栅格尺寸高

        // 2. 获取遥感影像经纬度范围，计算遥感影像像素分辨率
        //   获取原始影像的地理坐标范围
//        double latMin = ori_transform[3];
//        double latMax = ori_transform[3] + (yCount * ori_transform[1]);
//        double lonMin = ori_transform[0];
//        double lonMax = ori_transform[0] + (xCount * ori_transform[1]);

        double lonMin = ori_transform[0];
        double latMax = ori_transform[3];
        double lonMax = lonMin + (xCount * ori_transform[1]) + (yCount * ori_transform[2]);
        double latMin = latMax + (xCount * ori_transform[4]) + (yCount * ori_transform[5]);

        // 原始图像东西方向像素分辨率
        double src_w_e_pixel_resolution = (lonMax - lonMin) / xCount;
        // 原始图像南北方向像素分辨率
        double src_n_s_pixel_resolution = (latMax - latMin) / yCount;
        // 获取源图像crs 获取源影像的坐标参考
        String src_wkt = dataset.GetProjectionRef();
        SpatialReference spatialReference = new SpatialReference(src_wkt);
        String value = spatialReference.GetAttrValue("AUTHORITY", 0) + ":" + spatialReference.GetAttrValue("AUTHORITY", 1);
        CoordinateReferenceSystem crs = null;
        try {
            crs = CRS.decode(value);
        } catch (FactoryException e) {
            e.printStackTrace();
        }

        for (int zoom = minZoom; zoom <= maxZoom; zoom++) {
            // 3. 根据原始影像地理范围求解切片行列号
            int tileRowMax = getOSMTileYFromLatitude(latMin, zoom);
            int tileRowMin = getOSMTileYFromLatitude(latMax, zoom);
            int tileColMin = getOSMTileXFromLongitude(lonMin, zoom);
            int tileColMax = getOSMTileXFromLongitude(lonMax, zoom);

            for (int row = tileRowMin; row <= tileRowMax; row++) {
                for (int col = tileColMin; col <= tileColMax; col++) {

                    // 4. 求原始影像地理范围与指定缩放级别指定行列号的切片交集
                    double tempLatMin = tile2lat(row + 1, zoom);
                    double tempLatMax = tile2lat(row, zoom);
                    double tempLonMin = tile2lng(col, zoom);
                    double tempLonMax = tile2lng(col + 1, zoom);
                    ReferencedEnvelope tileBound = new ReferencedEnvelope(tempLonMin, tempLonMax, tempLatMin, tempLatMax, crs);
                    ReferencedEnvelope imageBound = new ReferencedEnvelope(lonMin, lonMax, latMin, latMax, crs);
                    ReferencedEnvelope intersect = tileBound.intersection(imageBound);

                    // 5. 求解当前切片的像素分辨率(默认切片大小为256*256)
                    // 切片东西方向像素分辨率
                    double dst_w_e_pixel_resolution = (tempLonMax - tempLonMin) / 256;
                    // 切片南北方向像素分辨率
                    double dst_n_s_pixel_resolution = (tempLatMax - tempLatMin) / 256;


                    // 6. 计算交集的像素信息
                    // 求切图范围和原始图像交集的起始点像素坐标
                    int offset_x = (int) ((intersect.getMinX() - lonMin) / src_w_e_pixel_resolution);
                    int offset_y = (int) Math.abs((intersect.getMaxY() - latMax) / src_n_s_pixel_resolution);

                    // 求在切图地理范围内的原始图像的像素大小
                    int block_xsize = (int) ((intersect.getMaxX() - intersect.getMinX()) / src_w_e_pixel_resolution);
                    int block_ysize = (int) ((intersect.getMaxY() - intersect.getMinY()) / src_n_s_pixel_resolution);

                    // 求原始图像在切片内的像素大小
                    int image_Xbuf = (int) Math.ceil((intersect.getMaxX() - intersect.getMinX()) / dst_w_e_pixel_resolution);
                    int image_Ybuf = (int) Math.ceil(Math.abs((intersect.getMaxY() - intersect.getMinY()) / dst_n_s_pixel_resolution));

                    // 求原始图像在切片中的偏移坐标
                    int imageOffsetX = (int) ((intersect.getMinX() - tempLonMin) / dst_w_e_pixel_resolution);
                    int imageOffsetY = (int) Math.abs((intersect.getMaxY() - tempLatMax) / dst_n_s_pixel_resolution);
                    imageOffsetX = imageOffsetX > 0 ? imageOffsetX : 0;
                    imageOffsetY = imageOffsetY > 0 ? imageOffsetY : 0;

                    // 7. 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。
                    // 推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。
                    Band in_band1 = dataset.GetRasterBand(1);
                    Band in_band2 = dataset.GetRasterBand(2);
                    Band in_band3 = dataset.GetRasterBand(3);

                    int[] band1BuffData = new int[256 * 256 * gdalconst.GDT_Int32];
                    int[] band2BuffData = new int[256 * 256 * gdalconst.GDT_Int32];
                    int[] band3BuffData = new int[256 * 256 * gdalconst.GDT_Int32];

                    in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf,gdalconst.GDT_Int32, band1BuffData, 0, 0);
                    in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf,gdalconst.GDT_Int32, band2BuffData, 0, 0);
                    in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf, gdalconst.GDT_Int32, band3BuffData, 0, 0);

                    //  8. 将切片数据写入文件
                    // 使用gdal的MEM驱动在内存中创建一块区域存储图像数组
                    Driver memDriver = gdal.GetDriverByName("MEM");
                    Dataset msmDS = memDriver.Create("msmDS", 256, 256, 3);
                    Band dstBand1 = msmDS.GetRasterBand(1);
                    Band dstBand2 = msmDS.GetRasterBand(2);
                    Band dstBand3 = msmDS.GetRasterBand(3);

//                    // 设置alpha波段数据,实现背景透明
//                    Band alphaBand = msmDS.GetRasterBand(4);
//                    int[] alphaData = new int[256 * 256 * gdalconst.GDT_Int32];
//                    for (int index = 0; index < alphaData.length; index++) {
//                        if (band1BuffData[index] > 0) {
//                            alphaData[index] = 255;
//                        }
//                    }

                    // 写各个波段数据
                    dstBand1.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band1BuffData);
                    dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData);
                    dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData);
//                    alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData);


                    String path = outputDir + "/" + zoom + "/" + col + "/";
                    String pngPath = path+ row + ".png";
                    System.out.println(pngPath + "   切片完成");
                    // 使用PNG驱动将内存中的图像数组写入文件
                    Driver pngDriver = gdal.GetDriverByName("PNG");
                    File file = new File(path);
                    if (!file.exists()){
                        file.mkdirs();
                    }
                    Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS);

                    // 释放内存
                    msmDS.FlushCache();
                    msmDS.delete();
                    pngDs.delete();
                    file.delete();
                }
            }
        }

        dataset.delete();
    }


    /**
     * 瓦片坐标转经度
     *
     * @param x
     * @param z
     * @return
     */
    public static double tile2lng(double x, int z) {
        return x / Math.pow(2.0, z) * 360.0 - 180;
    }

    /**
     * 瓦片坐标转纬度
     *
     * @param y
     * @param z
     * @return
     */
    public static double tile2lat(double y, int z) {
        double n = Math.PI - (2.0 * Math.PI * y) / Math.pow(2.0, z);
        return Math.toDegrees(Math.atan(Math.sinh(n)));
    }

    /**
     * 根据经度获取切片规范下的行号
     *
     * @param lon
     * @param zoom
     * @return
     */
    public static int getOSMTileXFromLongitude(double lon, int zoom) {
        return (int) (Math.floor((lon + 180) / 360 * Math.pow(2, zoom)));
    }

    /**
     * 根据纬度获取切片规范下的列号
     *
     * @param lat
     * @param zoom
     * @return
     */
    public static int getOSMTileYFromLatitude(double lat, int zoom) {
        return (int) (Math.floor((1 - Math.log(Math.tan(lat * Math.PI / 180) + 1 / Math.cos(lat * Math.PI / 180)) / Math.PI) / 2 * Math.pow(2, zoom)));
    }
}
